%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 模拟退火算法
% 简介：模拟退火算法寻找最优路径
% 作者：Zhaojiang
% 日期：2023/10/10
% 企鹅：277746470
% 最短长度：102.6600
% 最短路径：1->23->31->16->19->10->25->22->2->27->18->15->12->11->34->4->13->
%          20->17->24->14->5->21->8->29->7->9->30->28->26->3->6->32->33->1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;clear;clc;
%% target fucntion
tsp=TSP(importdata('./city_data.mat'));
figure('Name','模拟退火算法');
tsp.draw_path;
disp 当前路径
tsp.print_path;
pause(2);
close all;
figure('Name','模拟退火算法','Position',[200 200 1280 480]);
%% 设定算法参数
rng(1234);
max_iteration = 10000;% 最大降温次数
temp_init = 100;% 初始温度
temp = temp_init;% 当前温度
cooling_rate = 0.99;% 设置冷却速率
exchange_thre = 0.5;% 选择交换产生新解的概率
temp_end = 1;
%% 算法数据初始化
city_num = tsp.city_num;% 城市数量
best_iteration = 0;
best_path = tsp.path;% 最佳路径
best_length = tsp.path_length;% 最佳长度
all_best_length = zeros(1,max_iteration);% 记录最佳函数值变化
%% 开始迭代求解
for iteration=1:max_iteration
    % 生成新解
    old_path = tsp.path;
    new_path = old_path;
    if rand < exchange_thre
        % 交换位置
        randposition = randperm(city_num-1)+1;
        position1 = randposition(1);
        position2 = randposition(2);
        new_path(position1) = old_path(position2);
        new_path(position2) = old_path(position1);
    else
        % 交换位置
        randposition = randperm(city_num-1)+1;
        randposition = sort(randposition(1:3));
        position1 = randposition(1);
        position2 = randposition(2);
        position3 = randposition(3);
        son_path1 = old_path(position1:position2);
        son_path2 = old_path(position2+1:position3);
        new_path(position1:position3) = [son_path2,son_path1];
    end
    % 新旧解差值（新-旧）
    old_length = tsp.path_length;
    tsp.update_path(new_path);
    new_length = tsp.path_length;
    delta = new_length - old_length;
    % 判断是否接受新解
    if (delta < 0) || (rand <= exp(- delta / temp))
        tsp.update_path(new_path);
    else
        tsp.update_path(old_path);
    end
    % 是否当前最优解
    new_length = tsp.path_length;
    if new_length < best_length
        best_length = new_length;
        best_iteration = iteration;
        best_path = tsp.path;
    end
    all_best_length(iteration) = best_length;
    temp = temp*cooling_rate;
    if mod(iteration,max_iteration/500)==0 || iteration ==1
        subplot(1,2,1);hold off;
        tsp.update_path(best_path);
        tsp.draw_path;
%         disp 当前路径
%         tsp.print_path;
        disp(['迭代次数：' num2str(iteration) '/' num2str(max_iteration)])
        disp(['最短长度:' num2str(best_length)])
        subplot(1,2,2);hold off;
        plot(all_best_length(1:iteration));
        drawnow;
%         pause(0.1);
%         break;
    end
end
%% 处理结果
clf;
subplot(1,2,1);hold on;
tsp.update_path(best_path);
tsp.draw_path;
% 打印结果
disp(['最优迭代次数：' num2str(best_iteration)])
disp(['最短长度:' num2str(best_length)])
disp 最佳路径:
tsp.print_path();
% 可视化结果
subplot(1,2,2);
plot(all_best_length)
save('./Best_Path_SA.mat','best_path','-double');